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^ , We report the production and benchmarking of several refinements of the power method that 
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we used an observation by Booth that has made possible the calculation of up to the lO*'^ eigenpair 
for simple test problems simulating the transport of neutrons in the steady state of a nuclear reactor. 
Here, we summarize our techniques and efforts to-date on determining mainly just the two largest 
or two smallest eigenpairs. To illustrate the eff'ectiveness of the techniques, we determined the two 
extremal eigenpairs of a cyclic matrix, the transfer matrix of the two-dimensional Ising model, and 
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Qj , the Hamiltonian matrix of the one-dimensional Hubbard model. 
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I. INTRODUCTION 

Computing eigenpairs of large matrices is a ubiquitous problem in computational physics. 
In this paper, we present several refinements of the basic power method that enable the 
efficient and accurate computation of multiple extremal eigenvalues of very large matrices. 
Ultimately, our objective is producing Monte Carlo versions of such methods for matrices 
whose orders are so large that even the eigenvectors cannot be stored in computer memory. 
For such problems, the computation of a basic vector quantity as the inner product is 
generally either very inefficient or impractical. It can be impractical, for instance, because 
the nature of Monte Carlo sampling means most components of these vectors are unknown. 
Here, we focus on the basic algorithms developed to date, noting they work well when used 
deterministically. Novel will be the illustration of how the power method can be expanded 
to compute several extremal eigenpairs simultaneously rather than just one at a time. While 
various versions of the power method often compute very well the dominant eigenvalue Ai, 
the one with largest absolute value, computing subdominant eigenvalues A2, A3, . . . has often 
proven much more difficult and is much less frequently attempted. 

The algorithms to be presented use some recent insights of Booth [l|, |2| that were de- 
veloped for Monte Carlo simulations of steady state neutron transport in nuclear reactors. 
Initially, he proposed a novel modification of the power method that has produced up to 
10 eigenpairs for simple test problems. Here, we present refinements of these insights and 
focus on determining just Ai and A2, plus their eigenvectors. The convergence of the power 
method is well known to slow as the ratio A2/A1, sometimes called the dominance ratio, 
approaches unity. As such, this ratio is an indicator of solution difficulty and acceptability. 
As we will illustrate, an advantage of computing two dominant eigenpairs simultaneously is 
often improved convergence to the first one. An advantage of the present techniques is the 
ease in getting both eigenfunctions along with their eigenvalues. 

It is important to note that various areas of science and engineering seek multiple eigen- 
pairs for reasons other than algorithmic gains. In nuclear engineering, a dominance ratio 
distinct from unity is an acceptance qualifier for various nuclear criticality safety assessments 
and nuclear reactor designs 3|. In statistical physics, a dominance ratio nearing unity, on 
the other hand, is often a condition sought. Near a continuous phase transition, A2 — > Ai, 
and A2/A1 controls the microscopic spatial correlations among physical degrees of freedom 
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nomena, phase transitions driven by zero-point motion at zero temperature [51]. Here, it is 
the two smallest eigenvalues of the Hamiltonian matrix describing the physical system that 
are of interest. The quantum critical phenomenon construct, while supplemented by a few 
exact solutions to some very simple problems, is largely phenomenological in part because 
of the inability to compute A2 for models of direct physical relevance. 

In the next section. Section II, we summarize some basic notions about the power method 
and our refined procedures. For simplicity, we will assume the two largest extremal eigen- 
pairs are sought. Also we restrict attention to systems with real eigenvectors and eigenvalues, 
but our methods can be applied to complex systems as well. In Section III, we apply these 
techniques to determination of a few eigenpairs of three problems. The first is the cyclic 
matrix that results from the discretization of the gradient operator on a circle, the second is 
the transfer matrix of the two-dimensional Ising moel, and the third is the Hamiltonian ma- 
trix of the one- dimensional Hubbard model. For the first and third problems, we determine 
the smallest two eigenpairs (ground-state and first excited state) instead of the largest ones 
to illustrate the flexibility of the techniques. The second and third problems counterpose 
in their computational challenges in a number of ways: The transfer matrix for the Ising 
model is non-symmetric, positive, and dense. Its eigenvalues are known analytically and 
all its matrix elements follow a single simple analytic expression. The Hamiltonian matrix 
for the Hubbard model is symmetric, indefinite, and sparse. Its two smallest eigenvalues 
are not known analytically, and its matrix elements, while easy to compute, lack a simple 
expression. In Section IV, the final section, we will discuss extensions of the techniques to 
broader classes of problems, including those involving continuous operators. 

II. METHODOLOGY 

We first summarize the power method, and then we discuss ways to refine it so convergence 
is to the two largest extremal eigenpairs simultaneously. We conclude this section with two 
refinements of the power method: one is necessary for the reduction of round-off error and 
the other improves the convergence rate to the dominant eigenpair while simultaneously 
calculating the second extremal eigenpair. 



A. Power Method Basics 



For some real- valued N x N matrix A, not necessarily symmetric, we will be concerned 
with the N eigenpairs (Aj,'?/'i) satisfying 



A^l^i = Xi^pi 



(1) 



In the simplest application of the power method 6|], an iteration is started with some ip, 
normalized in some convenient, but otherwise relatively arbitrary, manner and consists of 

iterating two steps 

(j) = Alp 

^ = 0/11011 
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Hence, the dominant eigenpair is simultaneously determined. For the norm of the vector (p 
whose components are (pi, a frequent choice is 

II 4> 11 = 11 l|oo= max|0i| 

This is the choice adopted here. 

Clearly, if IA2/A1I ^ 1 convergence of the iteration is slow. Often it can be improved 
by the replacement A —>■ A — cxI which shifts the value of each eigenvalue by a constant 
amount a but does not change the associated eigenvector. Besides potentially accelerating 
convergence, the shift also enables the determination of the smallest, instead of the largest, 
eigenpair. In particular, if A and all the Aj are real, no matter how a is chosen, either Ai — cr 
or Aat — 0" will be the converged eigenvalue. Most often, a is chosen to be independent of 
iteration step. In this case, for convergence to Ai, the optimal choice for cr is | (A2 + \n)', 
for convergence to Aat, the choice is | (Ai + Aat-i) 6|]. 



If the dominant eigenvalue is degenerate, for example, doubly degenerate with Ai = A2, 
or degenerate in magnitude, for example, doubly degenerate with |Ai| = IA2I, then the power 
method, as most iterative methods, cannot determine a unique eigenvector. As can be seen 
from Eq. ([3]), in these situations the iteration converges to 



A"V^ = A^ 



ai^i + sign I — I ^2^2 + 2^ "i I T- I Vi 



The eigenvalue estimators will converge to the correct values of Ai and A2 but the eigenvector 
estimate corresponding to the dominant eigenvalue will be some linear combination of ipi and 
iIj2- a similar situation will can occur for convergence to the first subdominant eigenvalue if 
for example IA2I = IA3I. In this case ipi can be determined but ip2 cannot. 

If a few dominant eigenpairs, say M, are desired, one of two approaches are tried. One 
approach is to use the power method to determine the dominant eigenpair, use deflation 
to project out this state out of the matrix, and then reuse the power method on the de- 
flated matrix. To determine several eigenpairs simultaneously, the power method can be 
generalized to 

$ = A^ 

where $ and \I' are M x N matrices whose columns are orthogonalized to each other. 
This orthogonality needs maintenance throughout the computation or else all M vectors, 
represented by the columns of the initial \1/, will converge to the same one. This algorithm 
is called the simultaneous iteration method [7|. 

B. Observation 

Booth's reflnement of the power method [l|, l2| uses the observation that for any eigenpair 
(A, ip) and for each non-zero component of the eigenvector, the eigenvalue equation Aip = Xip 
can be rewritten as 

E Aa/3lpl3 

and that similar equations can also be written for any number of groupings of components, 

E E^a/3V'/3 E E^a/3V'/3 E E ^a/3^/3 
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where the Ri are rules for different groupings. The groups can overlap. In addition, any two 
groupings, say 1 and 2, imply 

E ^° H 5Z ^"/3^/? = E ^" E E ^c^pi^p (6) 

The eigenvalue estimator (j4]) is a special case of what is often called a mixed estimator |8| 

(01^) 

In the present case, the component (pi of the vector (p is unity if z G -R; otherwise, it is zero. 
From A^ groupings of the components. Booth constructs A^ estimators for the A^*'* eigenvalue 
and forces them to become equal by adjusting certain parameters at each iteration step. 
Several ways to do this have been devised, and we will now sketch the most recent ones for 
obtaining two extremal eigenpairs simultaneously. 

For almost any starting point ip = J^iO^ii^i, the power method will converge to {Xi,ipi). 
The same would be true for almost any other two normalized, but not necessarily orthogonal, 
starting points ip' = J2i hi^i or ip" = J2i cai^i- We will in fact choose two such starting points 
and at each step apply A to them individually, but at each step we will adjust the relationship 
between them to prevent the collapse of their sum to the dominant eigenf unction. 

Formally, we start the iteration with ip = ip' + rjip" and assume that after a large number 
of steps just the two dominant eigenpairs remain significant. Then we have 

2 2 

A> = A" ^ a,V. = E («^ + Vh) Xy^ (7) 

To determine r], we define two groupings of the components of A^ip, Ri and i?2, and let Kj 
be the eigenvalue estimate for the j*^ grouping. Then from Eq. ^ we find that 
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(ai + r]bi) A^ ^ ^L" + («2 + vh) A2 E ^2,a 
(ai + rih) Ar^ E ^i>" + («2 + Vh) ^2'' E ^2,a 

(ai + rih) \1 E V'l.a + (^2 + Vh) A2 E V^2,a 
ag-R2 ag-R2 

{a^ + r]h) X^' E ^i,« + («2 + ^^2) X^' E ^2,a 

a£R2 06-R2 



If we require ki = K2, a quadratic equation for 77 results. If one solution of this equation 
is chosen to guide the iteration to ai + r^foi = 0, then ni = k,2 = A2. If the solution on the 
other hand is chosen to guide the iteration towards 02 + 7762 = 0, then ki = K2 = Ai. 



In practice, we find the coefficients of this quadratic equation in the following manner: 
Suppose at the n*^ step, ip' and ip" have iterated to ip' and ^", then at the (n+ 1)*'^ step we 
require 



E V-^ + r] E < ~ E < + r7 E < 

oG-Ri oG-Ri aGR2 a£R2 

which leads to q2rf + qifj + go = with 
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50 = E < E E^«/3^^ - E < E E^"/3^/^ 
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The strategy is to apply A repeatedly until two real solutions for t] exist. One solution will 
then guide the iteration to (Ai,'?/'i); the other, to (A2,V^2)- Typically, this procedure would 
be used only if the simultaneous convergence to two pairs is desired or if the convergence to 
just the second eigenpair is desired. In some cases, however, accelerated convergence to the 
first pair occurs. 

C. First Refinement 

For simplicity, we focus on the determination of the second largest eigenpair [2] and note 
that one additional improvement is necessary for finite precision computers. As both ip' 
and ip" are converging to the first eigenfunction, only their sum, for proper choices of rj, is 
converging to the second one. Eventually, when rj is the root, say 772, guiding ip" toward the 
second eigenvector tp2, the determination of tp2 is limited by the accuracy of the sum of ip' 
and 772V'"- To mitigate this situation, we modify the iteration by making the replacements 
ip' <— tp' and ip" <— tp" + 7]2'ip' before moving to the {n + 1)*^ step, and then in the {n + 1)*'* 
step we find the new 77 from the quadratic equation and subtract from it the ri2 from the n*'' 
step. Formally, this is equivalent to rewriting the coefficients of the ipi in Eq. ([7]) as 

ai + biT] = {tti + bir]2) + bi{r] - r]2) (10) 



making the replacements 

Oj ^ (Oi + bi1]2) 

r] ^r]-r]2 

and then in the next iteration solving the quadratic equation for the shifted 77. 

What does this procedure accomplish? We note that near convergence, when only ipi 
and ■?/'2 remain significant, the current best estimate of ipi is contaminated with 1^2 and vice 
versa. Denoting these estimates by ipi + e-?/'2 and 1^2 + S'^Pi and introducing the adjustable 
parameter 77, we can write another estimate of ip2{v) as 

^2(^) = {i'2 + #1) + vii'i + (^^2) 

and so with the application of A to move to the (n + 1)*^ step, the new estimate of the 
second eigenfunction becomes 

i^2new = Aij2{r]) = {\2iJ2 + SXi^pi) + r]{\iipi + eA2V'2) 

If at this step 772 is the choice that guides to ■ip2, then we define 77" by 

r] = 1]" + 772 

so that 

1p2newiv') = (-^2^2 + ^Ai^i) + {t]" + 7/2)(Ai^i + eA2^2) 

= (A2^/'2 + (5AiV'i) +772(AiV^i + eA2V^2) + V(-^i^i + eA2^2) 

We observe that Xiipi+eX2ip2 is this step's power iteration estimate for the first eigenfunction 
so that the second eigenfunction, the term in brackets, is essentially being corrected by an 
attempt to remove the remaining contamination SXiipi from the first eigenfunction. If we 
define the new eigenfunction iterates as 

'ip2new = A2V'2 + SXi^pi 
'iplnew = XiiJi + eA2t/'2 

then 



Thus, the effect of Eq. ( ITOj) is promoting the convergence of ip' to the first eigenfunction in the 
normal power method way whereas tp" (the second eigenfunction estimate) is being corrected 
at each step by adding (removing) a httle of the first eigenfunction estimate. Convergence 
is reached when ^72 — * 0; that is, when the second eigenfunction needs no correction from 
the first eigenfunction. 

The above analysis leads to a simple numerical algorithm. The basic steps are 

Step 1: Initialize 

1. Set convergence parameter e to a small value. 

2. Choose initial estimates ip' ~ ipi and ij)" ~ V's 

Step 2: Reset 

1. Normalize V^' ^ V'VII^'ll and ^" ^ V /WW 
Step 3: Execute power step 

1. Apply A to ip' and ip" and solve resulting quadratic balance condition (Eq. ([8]) ) 

2. If the roots are real, assign the roots tji and 772 to correspond to the largest and 
smallest (in magnitude) eigenvalue estimates respectively and then update via 

ijj' ^ AiIj' 

^" ^ Alp" + 7]2A'iIj' 

else update via 

ijj' ^ Alp' 
^" ^ A^" 

Step 4: Test for convergence 

1. If |?72| > e, go to Step 2 

Step 5: Terminate 

Eigenvalue estimates can be made by placing ip' and ip" in Eq- (E]) for the same or different 
Ri. These Ri can be the same or different from the two used to compute the g^. When the 
roots are complex, an alternative to Step 3.2 is the use of complex arithmetic. If it is used, 
then the ip' and ip" estimates in Step 3.2 are updated with complex r/j. 
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D. Second Refinement 

Because Eq. ([7]) will be almost true for large n, it yields a way to estimate the ipi. If Cj 
are normalizing constants and 772 is the root that gives A2, then from Eq. i^^ the {n + 1)*'* 
guess at ip2 is 

If rji ^ —02/62 is the root that gives Ai, then from Eq. ([7]) the {n + 1)*'' guess at ipi is 

These two estimates suggest using 

as the next iteration guess. Next, we insert 

into Eq. (jS]) and solve it for the two r] roots. Now we take for new estimates (with the Q 
providing normalization) 



■ (n+2) 



^P2^^Pr"=C2[A^P^-^'^\,=,, (11) 

V;i^V^i"+^)=Ci[AV^("+^)],=,, (12) 

We note that Eq. (Ilip is the same adjustment to the second eigenfunction estimate as in the 
first refinement, which uses the best estimate of tp2- Equation flT^ uses the best estimate of 
ipi instead of the power iterated estimate used in first refinement. Empirically, this second 
refinement simultaneously produces estimates of ipi converging as A3/A1 and estimates of 
iIj2 converging as A3/A2. In the appendix we demonstrate these rates of convergence for 
non-degenerate states. 

Incorporating this refinement requires only replacing Step 3 of the algorithm for the 
previous with 

Step 3: Execute power step 

1. Apply A to ip' and iJj" and solve resulting quadratic balance condition (Eq. ([8])) 

10 



2. If the roots are real, assigning the roots rji and ri2 to correspond to the largest 
and smallest (in magnitude) eigenvalue estimates respectively and then update 
via 

^" ^ Alp" + r]2A'ilj' 



else update via 



ij' ^ AiIj' 

Step 4: Test for convergence 

1. If \rj2\ > e, go to Step 2 
Step 5: Terminate 

E. Practical Algorithm 

In an actual implementation of these algorithms, monitoring convergence by |?72| < e 
is not the only choice. The more common way would be monitoring successive estimates 
of the Aj plus monitoring the residuals || Aip' — Xiip' \\ and || Aip" — X2'ip" \\- We also 
note the following alternative: As ip' converges to ipi and ip" converges to 1^2, % and g2 
converge to zero. In short, multiple criteria exist, leading to cross checks. Some recycle 
already computed quantities and are consequently quite efficient. Here is an algorithm for 
the second refinement suitable for implementation: 

Step 1: Initialize 

1. Set convergence parameters eo and £2 to small values. 

2. Initialize iteration index n = 0, 

3. Choose initial estimates for ip' and ip" , 

4. Choose the rules -Ri and R2 for grouping of iterated vector components. 

Step 2: Reset 

11 



1. Normalize ip' ^ ip' /\\ip'\\ and ip" ^ ip^WW 
Step 3: Execute power step 

1. Apply A to Ip' and ip", 

2. Solve resulting quadratic balance condition (Eq. [8]), 

3. Estimate eigenvalues using either rule (region) Ri or R2, 

4. If the roots are real, assign the roots r]i and 172 to correspond to the largest and 
smallest (in magnitude) eigenvalue estimates. For example, we will have 

E j:Ao.ptP'f, + Vl E E^a/3^^ 

E E Aap'ip'. + r]2 T, T, AaM 

. _ a€Ri 13 ae-Ri /3 

A2 — 



E V'^ + ^2 E K 



then update via 



ip' <- Alp" + r]iA'ip' 
ip" ^ Alp" + r]2Aip' 



else update via 



iP' ^ AiP' 
iP" ^ AiP" 



Step 4: Test for convergence 



1. If either |go| > ^o oi' \(l2\ > ^2, increment the iteration index, n ^ n + 1 and go 
to Step 2. 

Step 5: Terminate. 

When the roots are complex, an alternative to Step 3.2 is the use of complex arithmetic. If it 
is used, then the eigenvector estimates in Step 3.2 are updated with complex rji. The choice 
of rules is quite flexible. A rule may use one vector component selected randomly, a small 

12 



number of components selected randomly, all odd or even components, the first or second 
half of the vector, etc. We note that g2 goes to zero faster than qo, and when g2 becomes very 
small, then the quadratic equation numerically reduces to qirj + go = which is solved to 
get ?72. Essentially ^2 = means that the dominant eigenpair is known to machine accuracy, 
so it cannot be improved on further iteration. 

III. APPLICATIONS 



A. Cyclic Matrix 

To illustrate the effectiveness of the first refinement, we applied it to the symmetric NxN 
matrix 
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whose eigenvalues for any N are 

7« 



2 — 2 cos kr, 



4 sin 



2 '^n 



^^ with n = 0, 1, 2, . . . , iV — 1. Physically, the matrix represents the A^ point 



where A;„ 

discretization of the second derivative defined on a circle. We note that for A^ odd, all but 
the minimal eigenstate {n = 0) are doubly degenerate, while for A^ even, all but the minimal 
{n = 0) and maximal (n = A^/2) ones are doubly degenerate. Accordingly for even A^, 

= 70 < 7l = 7Af-l < ■ ■ ■ < 7Af/2-l = 7Ar/2+l < lN/2 = 4 

Table 1 reports the results of a deterministic computation of the second smallest eigenpair 
for a sequence of even A^. To generate it, all the eigenvalues of the matrix were shifted 
by subtracting four times the identity matrix and then getting the two largest magnitude 
eigenvalues oi A — 41. For the shifted matrix Ai = 70 = —4 and A2 = 71 = — 4cos^(7r/A^). 
Ai is thus seen as being independent of A^ and is not reported. The A2's in Table 1 are 
4 plus the power method's computation of second largest magnitude eigenvalue oi A — 41. 
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The iteration was stopped when the absolute value of the maximum difference between any 
component of the eigenvector in successive iterations was less than 10"^'^. We see remarkable 
agreement between the values determined by the power method and the exact analytic value 
is obtained even for largest possible A^ on our desktop computer. We converged accurately 
to the second smallest eigenvalue even though it is approaching the smallest one as A^ is 
increased and is itself degenerate. By it being degenerate, our eigenvector estimate is a linear 
combination of two eigenstates that depends on the starting conditions for the iteration. It 
is not something we can benchmark but it does approximate well Aip = \2ip- For Ri and 
i?2, we used the first and second half of the vector components. For N up to about 1000, 
starting vectors whose components were set randomly worked well. For starting vectors at 
A^ > 1000 we used the eigenvectors found at N/2 injected into the higher dimension via 
a(2j) ^ 0.756(j) + 0.256(j + 1) and a(2j + 1) ^ 0.256(j) + 0.756(j + 1). The coefficients 
were chosen to adjust for the fact that a(2j) is 1 unit from b(j) and 3 units from b(j+l) and 
a(2j + l) is 1 unit from b(j + l) and 3 units from b(j). 

B. Two-Dimensional Ising Model 

The two-dimensional Ising model is one of the few two-dimensional models of a system of 
many interacting degrees of freedom that has an exact solution for its thermodynamic prop- 



Iq 



erties. This solution, first constructed by Onsager [9|], shows that in the thermodynamic 
limit the model has a phase transition between an magnetically ordered (ferromagnetic) 
state at low temperatures and a magnetically disordered state (paramagnetic) at high tem- 
peratures. Onsager succeeded in calculating many of the properties of the model exactly, 
including the temperature T^ at which the transition occurs. Key to his calculations was 
expressing the partition function of the model in terms of its transfer matrix [lOJ, finding 
the dominant eigenvalue of this matrix, and showing in the thermodynamic limit (letting 
the area of the model approach infinity) that this eigenvalue implies the onset of long-range 
ordering among the spin variables of the model. 

We will consider the model for finite area, that is, an m x n model defined with periodic 
boundary conditions in one direction and open boundary conditions in the other. Because 
of the one open boundary, the transfer matrix will thus be non-symmetric. In the absence 
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TABLE I: For the cyclic matrix, the exact and first refinement calculations for the sub-dominant 
eigenvalue A2, plus their difference. 



N 

100 

200 

400 

800 

1600 

3200 

6400 

12800 

25600 

51200 

102400 

204800 

409600 

819200 

1638400 

3276800 



Exact 
0.0039465433630297 
0.0009868793234571 
0.0002467350504105 
0.0000616847138537 
0.0000154212379171 
0.0000038553131951 
0.0000009638285310 
0.0000002409571473 
0.0000000602392878 
0.0000000150598221 
0.0000000037649555 
0.0000000009412389 
0.0000000002353098 
0.0000000000588274 
0.0000000000147069 
0.0000000000036766 



PMl 

0.0039465431649277 
0.0009868792721619 
0.0002467351362063 
0.0000616847980273 
0.0000154212887717 
0.0000038553389818 
0.0000009638405936 
0.0000002409625672 
0.0000000602416170 
0.0000000150608281 
0.0000000037654755 
0.0000000009414407 
0.0000000002353788 
0.0000000000589155 
0.0000000000148614 
0.0000000000036855 



Difference 
1.98E-10 
5.13E-11 
-8.58E-11 
-8.42E-11 
-5.09E-11 
-2.58E-11 
-1.21E-11 
-5.42E-12 
-2.33E-12 
-l.OlE-12 
-5.20E-13 
-2.02E-13 
-6.91E-14 
-8.82E-14 
-1.55E-13 
-8.88E-15 



of an applied magnetic field, the model's energy is 

m—l n m n 

i=l j=l i=l j=l 

Here, {i,j) are the coordinates of a lattice site. The Ising spin variable fiij on each site has 
the value of ±1, J > 0, and /ij,n+i = f^i,i- A column configuration of Ising spins will be 
denoted by 

j v/^lj' r''2j5 • • • ? H"m,j) 

and there are 2"* possible configurations for each column. 

The definition of the transfer matrix follows from the expression for the partition function 
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Z{m,n) = Eexp[-f3E{{fi})] 



E exp 

(Ti,...,cr„ 



-(3{E{Vi{a,) + V2{a,,a,+,)} 



where 



is the interaction energy of the j column and 



E L{ai, a2)L{a2, 0-3) ■■ ■ -L(a„_i, an)L{an, cxi) 

cri,...,(T„ 



m—l 

i=l 

th 



m 



V2 (o-j, aj+i) = -z^ 51 ^^hj^^ 



«,i+i 



i=l 



is the interaction energy between the j*^ and (j+1)*'^ columns, z/ = J/ksT, ks is Boltzmann's 
constant, and L{a, a') is the transfer matrix of order A^ = 2™ x 2" whose elements are 

(m— 1 \ / m \ 

z/ ^ /ifc/ifc+1 J exp ( i/ ^ /ifc/i^ J 

More succinctly, 

Z{m,n)=TTiL^) = Y.X] 

Onsager found analytic expressions for all the eigenvalues of the transfer matrix. Since 
we needed to form this matrix, we found it as convenient to compute them numerically. Spot 
checks produce excellent agreement between the the two approaches. In the thermodynamic 
limit, when T -^ T^, A2 — > Ai. Here, although we chose T = T^ and the order of or 
matrix became quite large, we were still reasonably far away for this critical point. Table HTl 
presents a comparison of the two largest eigenvalues of the transfer matrix as determined 
by our second refinement of the power method and those determined by the EISPACK 
eigensolver RG [11]. For cases except m = 11, we simply used as the solution the results 
after 100 iterations. For m = 11, we used 1000 iterations. The larger number of iterations 
was necessary to obtain the same level of accuracy. 

C. One Dimensional Hubbard Model 

The Hubbard Hamiltonian was originally proposed as a model for metallic ferromagnetism 



12l |. Most recently, its two-dimensional version has been the subject of intense scrutiny as a 
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TABLE II: Comparison of the two largest eigenvalues of the transfer matrix for the 2D Ising 
model as a function of the lattice edge size m computed by the EISPACK routine RG and the 
second refinement. N = 2"^ is the order of the matrixAlso shown is the fractional difference 
(Arg — ApM2)/A_RG between the different estimates. We choose m = n and T = Tc- 

PM2 FD 

3.41421355573626E+00 O.OOE+00 

1.41421355573626E+00 O.OOE+00 

7.46410158611907E+00 3.57E-16 

4.82842709270073E+00 -1.29E-15 

1.78770541980345E+01 -7.95E-16 

1.35518083939891E+01 3.93E-16 

4.41298558292434E+01 4.83E-16 

3.60398703210878E+01 5.91E-16 

1.10192319565854E+02 9.03E-16 

9.38962258961221E+01 -4.54E-16 

2.76599914093667E+02 -8.22E-16 

2.42266413140723E+02 O.OOE+00 

6.96269201662782E+02 1.47E-15 

6.21748520715909E+02 1.65E-15 

1.75565374661531E+03 1.30E-16 

1.59043428136461E+03 6.05E-12 

4.43180239838646E+03 -1.44E-15 

4.05958858259756E+03 1.79E-15 

1.11957434253463E+04 1.46E-15 

1.03466429299731E+04 8.79E-16 

2.82985308867954E+04 -3.60E-15 

2.63419326613632E+04 -3.87E-15 



m 


N 


RG 


1 


4 


3.41421355573626E+00 
1.41421355573626E+00 


2 


16 


7.46410158611908E+00 
4.82842709270073E+00 


3 


64 


1.78770541980345E+01 
1.35518083939891E+01 


4 


256 


4.41298558292434E+01 
3.60398703210879E+01 


5 


1024 


1.10192319565854E+02 
9.38962258961220E+01 


6 


4096 


2.76599914093667E+02 
2.42266413140723E+02 


7 


16384 


6.96269201662783E+02 
6.21748520715910E+02 


8 


65536 


1.75565374661531E+03 
1.59043428137424E+03 


9 


262144 


4.43180239838645E+03 
4.05958858259757E+03 


10 


1048576 


1.11957434253463E+04 
1.03466429299731E+04 


11 


4194304 


2.82985308867953E+04 
2.63419326613631E+04 
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possible model for electronic superconductivity. In one-dimension, a variant of it, called the 
Pariser-Parr-Popple Hamiltonian is frequently used to model conjugated cyclic molecules. 
Other variants model one-dimensional organic conductors. Because of the enormous amount 
of computer memory required by deterministic methods, precise specification of the ground 
state (zero temperature) properties of these models has often been hampered by techniques 
limited to relatively small system sizes. The memory requirements scale as 4^, where A^ is 
the number of lattice sites. 

The Hamiltonian operator for the Hubbard model is 

{i,j),o- i 

where the summation is over nearest-neighbor pairs of lattice sites i and j and electron 
spin (j; t and U are the hopping amplitude and repulsive Coulomb parameters; c^ „., Ci^„ and 
f>'i,a = H^Ci^a- are the creation, destruction, and number operators for an electron at site i 
with spin a. Usually a Fock basis is used to represent the Hamiltonian operator as a matrix 

hij = {i\H\i) 

where 

K) = l^i,T' ^2,T, • • • , «jv,t) I^i,i5 ^2,1, • • • , nN,i) 

with n°" = 0, 1 being the eigenvalues of the number operator. As a representative of an Her- 
mitian operator whose matrix elements are real, the resulting matrix is symmetric. Various 
symmetries are usually used to block diagonalize the matrix and then obtain the ground state 
for each block. We will consider only blocks that have a specific value of the z-component 
of the total electron spin. The size of the Hilbert space and hence the order of the matrix is 

N '"■ 



N,I{N-Ni)l 

where iV| is the number of up spin electrons and being and Ni is the number of down spin 
electrons. 

For a given lattice site i the maximum number of non-zero values of hij is 2zN where 
z is the number of nearest neighbors of the chosen lattice. Typically, z <^ N; hence, the 
matrix is very sparse. Here, we will consider the model in one-dimension where z = 2. In 
one-dimension the model has an exact solution. Obtaining the ground or first excited state 
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for these solution is not as straightforward as for the two previous test cases. We chose 
to obtain them numerically and compare the effectiveness of our second refinement of the 
power method to that of several standard eigenpair methods. Our emphasis is on how well 
degenerate states are captured. 

We will compute the two largest and two smallest eigenvalues of the sparse, potentially 
hugely dimensioned, matrix H representing the operator H in one dimension with periodic 
boundary conditions. In this case and if N^ = Ni and N^ is odd, the model satisfies the 



following version of the Perron- Frobenius Theorem [l3j: If a matrix is irreducible and all 
off-diagonal elements are non-positive, the state corresponding to smallest eigenvalue is real 
and non-degenerate. We will study the model on a 10 site lattice with U = 4 and t = 1. For 
this lattice size and filling half or less, the theorem applies to cases {N-i,Ni) =(1,1), (3,3), 
and (5,5). They are called closed shell cases and the result of our calculations for them 
are shown in Table IIIII Because of various symmetries, features of the smallest states are 
reflected in those of the largest ones. 



In Table IIIII three methods where used to get the eigenvalues. One method used was 

n 

the LAPACK routine DSYEV |14|. This double precision routine returns all the eigenvalues 
and eigenvectors of a symmetric matrix. At large orders computer memory became insuffi- 
cient for its use. Accordingly, we supplemented our results with those obtained by using the 



DNLASO double precision subroutine [15| which is a block Lanczos method with selective 
reorthogonalization [7]. The components of the starting vectors for the Lanczos iteration 
are selected randomly and uniformly on the interval (-0.5,0.5). The quality of the results is 
controlled by specifying the block size (the number of starting states), the number of signif- 
icant figures for the convergence of the eigenvalues, and the maximum number of iterations. 
We found a block size of 1 gave estimates for the second and third eigenvalues that became 
progressively poorer as A^ increased. This is reasonable. A size of 2 produced cases where 
the sub-dominant eigenvalue was consistently returned as the dominant. For a size of 6, 
convergence was very slow if at all. For a size of 8, memory soon became insufficient. A 
size of 4 was used for the data in the table. We found lack of convergence for many cases 
if the precision was requested to be larger than 8 decimal places. Typically, a few hundred 
iterations were needed, but the computation times were a few tens of seconds. Table UlI] 
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TABLE III: For closed shell cases, comparison of the eigenvalues of a 10 site ID Hubbard model 
computed by the eigenpair routine DSYEV, the block Lanczos routine DNLASO, and the sec- 
ond refinement. For the first two methods, the three largest and three smallest eigenvalues were 
computed to measure their consistency, effectiveness, and accuracy. 

iV^ iV| N DSYEV DNLASO PM2 

1 1 100 0.5657693716217906E+01 0.5657693716217901E+01 0.5657693716217914E+01 

0.5519554669107880E+01 0.5519554669107876E+01 0.5519554669107137E+01 

-0.3862202348191250E+01 -0.3862202348191248E+01 -0.3862202348191251E+01 

-0.3618033988749895E+01 -0.3267468797160054E+01 -0.3618033988603501E+01 

3 3 14400 0.1656339684606611E+02 0.1656339684376816E+02 0.1656339684606624E+02 

0.1617312172182284E+02 0.1617312172136987E+02 0.1617312172191405E+02 

-0.8262531385370846E+01 -0.8262531383972004E+01 -0.8262531385370927E+01 

-0.7599976793651736E+01 -0.7599976793264113E+01 -0.7599976793831864E+01 

5 5 63504 0.2583432263352126E+02 0.2583432263577081E+02 

0.2543485463377173E+02 0.2543485464252857E+02 
-0.5834322635176973E+01 -0.5834322635773042E+01 
-0.5434854632148166E+01 -0.5434830052960784E+01 

shows excellent agreement between all three methods. We note that the excited state was 
at least doubly degenerate. 



The next set of results are for electron fil 



Kramers degeneracies. Kramers's Theorem IGj says all energy levels of a system containing 



ings where the eigenstates are subjected to 



an odd number of electrons must at least be at least doubly degenerate provided there are no 
magnetic fields present to remove time-reversal symmetry. In Table IIVI we present several 
Kramers cases where {N-^, Ni) = (3,2), (4,3), and (5,4). For the standard software packages, 
we listed the three largest and three smallest eigenvalues they estimated to see the accuracy 
to which they determined the degenerate ground state. We see very good agreement between 
all three methods in estimating the eigenvalue of degenerate largest and smallest state. All 
three however lack the precision necessary to differentiate between a true degeneracy and a 
very near one. The power method in particular is less than adequate for this purpose. 
In Table IIVI the three largest and smallest eigenvalues are presented to provide extra 
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information about the degeneracies. Sixteen significant figures were printed to indicate 
how well degeneracies are captured. Basically we do not know how the precision of the 
eigenvalues other then it is no more than difference between eigenvalues that should exactly 
be degenerate. All the expected features of the lowest eigenstates with regard to degeneracies 
are exhibited. Because the model has particle-hole symmetry similar features also exists for 
its largest eigenstates. 



For other electron fillings a priori exact information about degeneracies is lacking. What 
is known is that when f/ = the ground state for most fillings is degenerate. These fillings 
typically are called open shell cases. When f/ 7^ 0, the degeneracies are typically lifted, the 
degree to which it is however depends on the closeness of the nearest unoccupied eigenstate. 
For small systems, the lifting might be minor, and for limited precision calculations this 
can make distinguishing degenerate and nearly degenerate states difficult. For the results 
in Table IVl this discussion applies to the cases {N^.N^) = (2,2) and (4,4) illustrated there. 
The qualitative character of the results are similar to those of the Kramers's cases. 



The results of the three tables indicate that this test case is nontrivial. We judge our 
preliminary results as indicating that the block Lanczos and the second refinement of the 
power method have comparable effectiveness. The version of the power method used here 
lacks the ability to return more than two eigenvalues in contrast to the block Lanczos used 
that should estimate well at least four. 

We comment that the Lanczos method is not a black box. To ensure that the result 
is the minimum as opposed to some excited state, the calculation usually needs to be run 
multiple times with different random number seeds or some other means to change the 
starting vectors, and then the results need to be studied to identified those to be regarded 



as estimates of the ground state 



0, 



18j. An error is usually estimated from the variance of 
the average of the ground state estimates. We restarted the Lanczos calculations for a few 
of the cases multiple times. As our interest was qualitative, we did not perform an error 
analysis but instead presented representative results. 
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TABLE IV: For Kramers degeneracies cases, comparison of the eigenvalues of a 10 site ID Hubbard 
model computed eigenpair routine DSYEV, the block Lanczos program DNLASO, and the second 
refinement of the power method. For the first two methods, the three largest and three smallest 
eigenvalues were computed to measure their consistency, effectiveness, and accuracy in computing 
degenerate eigenvalues. 

iV| iV| N DSYEV DNLASO PM2 

3 2 5400 0.1306499556833340E+02 0.1306499556833341E+02 0.1306499556833335E+02 

0. 1306499556833336E+02 0. 1306499554321960E+02 0. 1306499556833335E+02 

0.1282579739183819E+02 0.1282579739163641E+02 

-0.7511951740365890E+01 -0.7511951740281242E+01 -0.7511951740365513E+01 

-0.7511951740365851E+01 -0.7511951731564561E+01 -0.7511951740365509E+01 

-0.7249884543021683E+01 -0.7249884543021723E+01 

4 3 25200 0.1816344283994604E+02 0.1816344283994595E+02 0.1816344283994610E+02 

0. 1816344283994604E+02 0. 1816344283994549E+02 0. 1816344283994610E+02 

0.1771746494384758E+02 0.1771746494296794E+02 

-0.8030089029893539E+01 -0.8030089029893475E+01 -0.8030089030622399E+01 

-0.8030089029893492E+01 -0.8030089029485101E+01 -0.8030089030532327E+01 

-0.7521441552342070E+01 -0.7521441551092671E+01 

5 4 52920 0.2285321122055267E+02 0.2285321221296537E+02 

0.2285321121726444E+02 0.2285321166226352E+02 
0.2221382316719896E+02 

-0.6853211221744196E+01 -0.6853211215825024E+01 
-0.6853211221310334E+01 -0.6853211214979680E+01 
-0.6213823170572150E+01 

The power method is also not a black box. It used the same sparse matrix as the block 
Lanczos method. From it we only show the two largest and two smallest as our double 
precision code did not incorporate a procedure to allow the determination of the third 
eigenpair. The components of one starting state were selected uniformly and randomly over 
(0,1); the other from (-0.5,0.5). We defined our regions in the following manner. First we 
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TABLE V: For open shell cases, comparison of the eigenvalues of a 10 site ID Hubbard model com- 
puted by the LAPACK routine DSYEV, the Netlib program DNLASO, and the second refinement 
of the power method. For the first two methods, the three largest and three smallest eigenvalues 
were computed to measure their consistency, effectiveness, and accuracy in computing degenerate 
eigenvalues. 

iV| iV| N DSYEV DNLASO PM2 

2 2 2025 0.1121466372028744E+02 0.1121466372028743E+02 0.1121466372028747E+02 

0.1096186919469933E+02 0.1096186919469927E+02 0.1096186919053599E+02 
0.1096186919469928E+02 0.1096186919469929E+02 

-0.6601239688910290E+01 -0.6431629846631373E+01 -0.6601239688910274E+01 

-0.6431629846631359E+01 -0.6424903541072491E+01 -0.6431629865616197E+01 
-0.6431629846631350E+01 -0.5854101965765123E+01 

4 4 44100 0.2143485463460406E+02 0.2143485463565059E+02 

0.2106806509093548E+02 0.2106806509410040E+02 
0.2106806508923771E+02 

-0.7647179205940428E+01 -0.7647179208191599E+01 

-0.7538791441444121E+01 -0.7538797509518616E+01 
-0.7538791440796984E+01 

performed a random permutation of the vector components. Region 1 was the first N/2 + 1 
of these permuted components; region 2, the last N/2 + 1. Initially, we used cq = 10^^^ 
and £2 = 10~^°. Multiple starts with changing random number seeds were used to ensure 
consistency. We first computed the largest two eigenvalues and used the largest value, 
truncated to two significant figures, as the shift to compute the two smallest. Cases where 
the lowest or the largest had a near degeneracy converged only so far. Often with |go| < ^o 
being satisfied while the €2 being too small for g2 to converge. Instead of adjusting these 
stopping criteria, we found that simply stopping the iteration after some fixed number of 
iterations and the choosing the result by locating when the residual WAip" — X2ip"\\ ceased 
decreasing was very effective. 
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IV. CONCLUDING REMARKS 

We presented two refinements of the power method that enable the simultaneous de- 
termination of two extremal eigenpairs of a matrix. We illustrated their effectiveness by 
benchmarking them on three quite distinct but physically challenging problems. For the 
cyclic matrix, we exactly knew the eigenvalues and their degeneracies. We showed we could 
determine the two smallest extremal eigenvalues to nearly machine precision. For the transfer 
matrix of the two-dimensional Ising model, we knew the exact values. The two-dimensional 
Hubbard model was more challenging: As a function of the electron filling various degen- 
eracies exist. Here, we choose to compare the effectiveness of the second refinement with 
two standard numerical determinations of the ground and first excited state, illustrating 
the limitations of all three methods especially when the dominant state is degenerate or 
very nearly so. In general, the second refinement appears as effective as a readily available 
implementation of the block Lanczos method with selective reorthogonalization. 

All our test cases involved real matrices, but preliminary testing indicates that the tech- 
niques presented also work well for complex matrices and non-symmetric matrices with 
complex eigenvalues. The techniques presented are easily adapted to the determination of 
just the dominant or just the sub-dominant eigenpair. Convergence to the dominant one is 
accelerated as it is controlled by A3/A1 is instead of A2/A1. Generalizations to more than 
two extremal pairs are possible. Preliminary testing for up to four have been promising. 

One advantage of our refinements is that they maintain the simplicity of the basic power 
method. Another is their adaptability to Monte Carlo implementations. In a number of 
fields of physics and chemistry, the power method is the core of the Monte Carlo meth- 
ods for determining the ground state of models whose complexity grow exponentially with 
physical size. Here, the ground state energy is estimated from samplings of the ground 
state wavefunction. Such samplings may involve only a small fraction of all possible compo- 
nents of the state, and a mixed estimator [S] for the energy is often used. Here we pushed 
the use of such estimators a step further with novel consequences. In two other papers, 
we will describe Monte Carlo implementations of our refinements and their applications to 
the transfer matrix of the two-dimensional Ising model and to Hamiltonian matrix of the 
two-dimensional Hubbard model [l9(]. 

The work was supported by the U. S. Department of Energy under the LANL/LDRD 
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APPENDIX 



We now show that if the eigenstates are non-degenerate, our second refinement of the 
power method simuhaneously produces estimates of ipi converging as A3/A1 and estimates 
of ■?/'2 converging as A3/A2. Suppose that the estimates of ipi and 1^2 are very good with only 
a small mixtures of other components; that is, there is some small v such that 

ipi^ {1 + dv)^i + ev^2 + fvi^s (A.l) 

^2 ~ av^i + (1 + bv)ij2 + cv^Js (A. 2) 

Hence 

■0 = {avtpi + (1 + bv)^lj2 + cv^ps) + x((l + dv)tlJi + evil)2 + fvip'i) (A. 3) 

Define the total i*'* eigenfunction component in region j as 

We now apply the balance condition for equal A's. For region Ri the eigenvalue estimate is 

NiiXi{av + x{l + dv)) + A"2iA2((l + hv) + xev) + N^iX^^cv + xfv) 
Nn{av + x{l + dv)) + N2i{{l + bv) + xev) + N^iicv + xfv) 

while for region R2 it is 

A^i2Ai(a^ + x{l + dv)) + iV22A2((l + bv) + xev) + N32X3{cv + xfv) 
Nuiav + x{l + dv)) + iV22((l + bv) + xev) + N32{cv + xfv) 

Set the above two eigenvalue estimates equal and cross multiply to clear the denominators 
and obtain a quadratic equation in x. We now collect powers of x. Terms involving x° are: 

Li = -a\iNi2N2iV + a\2Ni2N2iV + a\iNiiN22V - a\2NiiN22V (A.5) 

L2 = -c\2N22NsiV + c\3N22N3iV + c\2N2lNs2V - c\sN2lN32V 

L3 = v^{-abXiNi2N2i + a6A2A^i2A^2i + a6AiA^iiA'22) 

L4 = -CX2N22NS1V + CXSN22N31V + CX2N21NS2V - cA3A^2lA^32'y 

L5 = v'^{-bcX2N22Nsi + bcX3N22Nsi + acXiNnN32) 



Lq = v\-acXsNnNs2 + bcX2N2iN32 - bcXsN2iNi 



32 j 
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Terms involving x^ are 



Lr = -f\2N22NsiV + f\3N22N3lV + f\2N2lN32V - f\3N2lN32V (A.6) 

Lg = v'^{-ae\iNuN2i + aeA2iVi2iV2i) 

Lg = v^{aeXiNnN22 - aeX2NnN22 - afXiNuN^i) 

Lio = v^afXsNuNsi - ceAsATasATgi - fo/AsATasATgi) 

Ln = v\ceX3N22N3i + bfX3N22N3i + afXiNuN32) 

Ll2 = V^i-afX3NnN32 + CeX2N2lN32 + &/A2iV2lA^32) 

Ll3 = v\-CeX3N2iN32 - bfX3N2iN32) 

Li4 = (1 + dv){-XiNi2N2i + X2N12N21 + XiNuN22 - X2NUN22) 

Li5 = v{l + dv){-bXiNi2N2i + bX2Ni2N2i + bXiNnN22-bX2NnN22) 

Lie = v{l + dv){-cXiNi2N3i + cX3Ni2N3i + cXiNnN32-cX3NnN32) 



Terms involving x'^ are 



Ln = v^ef{-X2N22N3i + X3N22N31 + X2N21N32 - X3N21N32) (A.7) 

Lis = ev{l + dv){-XiNi2N2i + A2Ari2A^2i + AiA^iiA/'22 - A2A^iiA^22) 
Li9 = /t;(l + dv){-XiNi2N3i + A3Ari2Ar3i + AiiViiAr32 - A3ArnAr32) 

Finally we can write 

go = L1 + L2 + L3 + U + U + Le (A.8) 

qi = L7 + Lg + Lg + Lio + Ln + L12 + L13 + Lu + L15 + Lie (A.9) 

gs = L17 + L18 + L19 (A. 10) 

and note that we seek the smallest magnitude root of 

g2a:^ + gia; + go = (A.ll) 

Note that x will be very small when f ?« 0, so that the quadratic part can be ignored, leading 
to 

a: = -go/gi (A. 12) 

Also the terms involving f^ can be ignored compared to terms involving v. Thus 
go = Li + L2 
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= v{-a\iNi2N2i + aA2A^i2A^2i + a\iNiiN22 - a\2NnN22 
-CX2N22N31 + CX3N22N31 + CX2N21N32 - CX3N21N32) 

Now consider qi. Every term except L14 has at least v^ in it, which will be small compared 
to the f° term in L14. Thus using the f° term in L14, 

qi = {-X1N12N21 + A2A^i2A^2i + XiNnN22 - A2iViiiV22) 
so from Eq. (1A.12I) we have that 

[-X2N22N3I + X3N22N;^ + X2N21N32 - A3A^2lA^32) 



X ^ — dv — cv~ 

{-X1N12N21 + X2NUN21 + XiNnN22 - X2NnN22) 

In what follows, ip is replaced by ip2 because the root corresponding to ip2 has been selected. 

A'ijj2 = Xi{av + x{l + dv))il>i + A2((l + hv) + xev)il>2 + X'i{cv + xfv)il>'i 
Both X and v are small, so the product xv is ignored yielding 

A'4)2 = Xi{av + x)'iIji + A2((l + bv) + xev)'4>2 + A3(cf + xfv)^3 
Substituting for x, we rewrite this equation as 

-Cv{-X2N22Nu + X3N22N31 + X2N21N32 - X3N21N32) 

"^^ ' {-XiNi2N2i + X2NuN2i + XiNuN22-X2NnN22) 
+A2(1 + bv)^2 + X-i{cv)il)3 

Dividing by A2(l + few) and keeping terms to order v yields 

]_ ^ -Cv{-N22N-il + (A3/A2)Ar22A^31 + A^2lA^32 - (A3/A2) Ar2lA^32) , 

A2 "^^ (-iVi2Ar2i + (A2/Ai)A^i2iV2i + ArnAr22 - (A2/Ai)iVniV22) "^^ 

+^2 + -^{cv)il)3 

We note that the ip^ component has dropped by the ratio A3/A2 after the iteration. Note 
that the tpi component in the above is proportional to the ip^ component (cv) in the estimate 
of 1^2 at the beginning of the iteration. Therefore, both the ipi and ips components drop out 
of the iIj2 estimate as A3/A2. 

Now we will look at the estimated ipi component. We note that dividing the iterate by 
a constant does not affect the eigenvalue estimates so we divide ( ]A.3|) by x and then label 
the estimate as ipi because the root corresponding to tpi will be chosen. 

ipi = {avipi + (1 + hv)il)2 + cvipz)- + ((1 + dv)'4>i + evi>2 + fvip^) 

Jb 
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The eigenvalue balance equation is the same, but instead of Eq. lA.lll we have 

1 1 

When estimating ipi, note that 1/x will be very small when f ?^ 0, so that the inverse 
quadratic part of this equation can be ignored, leading to 

1 



X 



-q2/qi (A. 13) 



After all f ^ terms in g2 are dropped, flA.lOp and flA.SP become 



q2 = ev{-\iNi2N21 + A2A^12A^21 + \lNuN22 - \2NuN22) + 
fv{-XiNi2N3l + X3Ni2N3l + XlNuN32 - X3NUN32) 

and after all v and f ^ terms in gi are dropped, flA.lOp and flA.7l) become 



gi = i-XiNi2N2i + X2N12N21 + XiNnN22 - X2NnN22) 
From flATTSJ) 



1 ^ _^^ _ fv{-XiNi2N3i + X3N12N31 + XiNnN32 - AsATnATsa) 

X ^^ {-X1N12N21 + A2A^12A'21 + AiA^iiA'22 - A2A^llA'22) 

We now note that 

^Vi = {avXi'il)i + (1 + bv)X2ip2 + cvX^ips)- + ((1 + dv)Xiipi + eT;A2'i/'2 + f^X^ips) 

X 

If the small terms associated with f ^ and v- are dropped, then dividing the equation by Ai 
yields (to first order in v) 

—A^pl = -^i)2- + (^1 + ew-^^2 + fv^^s) 
Ai Ai X Ai Ai 

Substituting for -, 

, JX2/Xi){-XiNi2N3i + X3Ni2N3i + XiNuN32 - A3ArnAr32) 

(-AiA^12A'21 + A2Ari2Ar2i + AiA^iiA^22 - A2A^llA^22) 

We see that the ip3 component is decreasing as A3/A1 and that the ■ip2 component is propor- 
tional to the tps component {fv) at the beginning of the iteration. Thus the 1^2 component 
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also should be falling as A3/A1. Thus the procedure is converging to the first eigenfunction 
at the accelerated rate A3/A1. 
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